Chameleon f(R) gravity in the virialized cluster 
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Current constraints on f(R) gravity from the large-scale structure are at the verge of penetrating 
into a region where the modified forces become nonlinearly suppressed. For a consistent treatment 
of observables at these scales, we study cluster quantities produced in chameleon and linearized Hu- 
Sawicki f(R) gravity dark matter iV-body simulations. We find that the standard Navarro-Frenk- 
White halo density profile and the radial power law for the pseudo phase-space density provide 
equally good fits for f(R) clusters as they do in the Newtonian scenario. We give qualitative 
arguments for why this should be the case. For practical applications, we derive analytic relations, 
e.g., for the f(R) scalar field, the gravitational potential, and the velocity dispersion as seen within 
the virialized clusters. These functions are based on three degrees of freedom fitted to simulations, 
i.e., the characteristic density, scale, and velocity dispersion. We further analyze predictions for 
these fitting parameters from the gravitational collapse and the Jeans equation, which are found to 
agree well with the simulations. Our analytic results can be used to consistently constrain chameleon 
f(R) gravity with future observations on virialized cluster scales without the necessity of running a 
large number of simulations. 



I. INTRODUCTION 

Modifications of gravity can serve as an alternative ex- 
planation to the dark energy paradigm for the late-time 
accelerated expansion of our Universe. Here, we special- 
ize to f(R) gravity, where the Einstein- Hilbert action is 
supplemented with a free nonlinear function f(R) of the 
Ricci scalar R [lj]. It has been shown that such models 
can reproduce the cosmic acceleration without invoking 
dark energy 0tJ|- However, they also produce a stronger 
gravitational coupling and enhance the growth of struc- 
ture. f(R) gravity is formally equivalent to a scalar- 
tensor theory where the additional degree of freedom is 
described by the scalaron field Jr = df/dR [a, Q- We 
parametrize our models by the background value of the 
scalaron field today, |/ro|. The fa field is massive, and 
below its Compton wavelength, it enhances gravitational 
forces by a factor of 4/3. Due to the density dependence 
of the scalaron's mass, however, viable f(R) gravity mod- 
els experience a mechanism dubbed the chameleon ef- 
fect 0-Q , which returns gravitational forces to the stan- 
dard relations in high-density regions, making them com- 
patible with solar-system tests [l(| at r < 20 AU. 

The enhanced gravitational coupling can be utilized to 
place constraints on the f(R) modification. The transi- 
tion required to interpolate between the low curvature 
of the large-scale structure and the high curvature of 
the galactic halo sets the currently strongest bound on 
the background field, \f R0 \ < |*| - (1CT 6 - 1CT 5 ) [3, 
i.e., the typical depth of cosmological potential wells. A 
bound of the same order is obtained from galaxies serving 
as strong gravitational lenses [1J| at r ~ (1 — 10) kpc. In- 
dependently, strong constraints can also be inferred from 
the large-scale structure (r > 1 Mpc). An upper bound 
°f l/fiol ^5 10~ 3 , for instance, can be obtained from the 



cluster density profiles constrained by weak lensing mea- 
surements [i2|. The currently strongest constraints on 
f(R) gravity models from the large-scale structure are 
inferred from the analysis of the abundance of clusters, 
yielding a constraint of \f R0 \ < 10" 4 [Hill. 

It is important to note that these cluster-scale con- 
straints have been derived by relying on a linearized ap- 
proach of the f(R) modifications, i.e., assuming a linear 
relation between the curvature fluctuation SR and the 
field fluctuation 5 fa that is correctly described by the 
background Compton wavelength parameter. This ap- 
proach, however, breaks down for |/i?o| < 10~ 5 , where 
cluster scales are affected by the chameleon mechanism. 
It is therefore important for comparison to future mea- 
surements to describe the observable quantities encom- 
passing the chameleon effect. 

Dark matter TV-body simulations of f(R) gravity pro- 
vide a great laboratory for the study of the chameleon 
mechanism, and many efforts have been made in perform 
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ing such simulations. For example, Oyaizu et al 
performed A-body simulations of the Hu-Sawicki 
f(R) gravity model for the first time using a particle 
mesh code. Later Zhao et al. [17| and Li et al. [18[ sim- 
ulated the same model using an adaptive particle-mesh 
code and significantly improved the resolution. 

In this paper, we aim at finding simple analytic and 
semi-analytic descriptions for cluster characteristics pro- 
duced in f(R) A-body simulations in both the linearized 
and chameleon scenarios. The relations we find here in- 
corporate the chameleon mechanism and can be used to 
assist in the consistent comparison of f(R) gravity to ob- 
servations. The outline of the paper is as follows. In ?JTT1 
we review f(R) gravity with a particular focus on the 
Hu-Sawicki model. In £ )III[ we provide (semi-)analytic 
relations for the scalar field, the gravitational potential, 



and the velocity dispersion at the virialized scales of clus- 
ters produced in linearized and chameleon f(R) gravity. 
Thereby, we start from the assumption of a Navarro- 
Frenk- White (NFW) profile for the cluster density and 
a power- law pseudo-phase-space density (PPSD). J jIVI is 
devoted to the comparison of these relations to the out- 
put of f(R) gravity TV-body simulations of gravitation- 
ally interacting cold dark matter particles. The fit of the 
cluster quantities is done using three degrees of freedom, 
i.e., the characteristic amplitude and scale of the NFW 
profile, and the velocity dispersion at the characteristic 
scale, we compare the values of these fitting parameters 
to predictions from scaling relations based on the spher- 
ical collapse and an estimation of the amplitude of the 
velocity dispersion employing the Jeans equation. We 
discuss our results in SjVj In the appendix, we give fur- 
ther details on the radial dependence of the PPSD profile 
used for the derivation of the velocity dispersion based 
on the expectations from the self-similar infall of a colli- 
sional gas in the context of enhanced gravitational forces 
such as present in f(R) gravity. We also motivate the 
applicability of the NFW and PPSD profile for clusters 
produced with modified gravitational forces. 



II. f(R) GRAVITY 



where 



9^ = (1 + fB.)g^, 
3 1 



df R J - 2k 2 (l + f R y 

m s TO 

VW = 2«»(1 +/*)'■ 
Integration of Eq. (J5)) gives the scalar field 

= ij|ln(l + /«) + o , 



(4) 
(5) 

(6) 

(7) 

(8) 



where we set 4>o = 0- Variation of the action with respect 
to <f> yields 



U<t> = -at + V (cf>) = y e ' ff (0), 



(9) 



where a = d\nA/d(j) and V e g is an effective potential 
governing the dynamics of <j>. Note that T = A(4>) 4 T. 
In the quasistatic limit, we neglect time derivatives in 
Eq. © and we obtain the scalar field equation of interest 

here. 



V*<t> = aA{4>fp m + V'(<j>), 



(10) 



where we assumed matter dominance and use physical 
coordinates. 



In f(R) gravity, the Einstein-Hilbert action is supple- 
mented by a free nonlinear function of the Ricci scalar 

R, 



S=—i <Fx^f—g [R + f(R)} + S m (^ m ; g^) . (1) 



Here, k 2 = 8-7T G, S m is the matter action with matter 
fields ip m , and we have set c = 1. Variation with respect 
to the metric g^ v yields the modified Einstein equations 
for metric f{R) gravity, 

G^+f R R^-({- - Df R ) g^-V^„f R = k 2 T^, (2) 



where the connection is of Levi-Civita type and f R = 
df/dR is the additional scalar degree of freedom of the 
model, characterizing the force modifications. 

By a conformal transformation of the metric, the Jor- 
dan frame action Eq. (fTJ can be recast in the Einstein 
frame. 



-S B1 [tp m ;A 2 (4>)g f _ ll , 



R 



lx, 



— --frw^-vw 



(3) 



Hu-Sawicki model 

We specialize our considerations to the functional form 
of f(R) proposed by Hu & Sawicki flQI ] . 



f(R) = -» 



ci (R/m 2 ) n 



(11) 



c 2 (R/m 2 ) n + 1 ' 

where m 2 = k 2 p m /3 and overbars refer to background 
quantities. The free parameters of the model c\ , c 2 , and n 
can be chosen to reproduce the ACDM expansion history 
and satisfy solar-system tests [lj| through the chameleon 

mechanism [7H9|]. In the high-curvature regime, c 2 R ^ 



2 , Eq. (|TT]) simplifies to 



JK ' c 2 n R n ' 



(12) 



where Rq denotes the background curvature today, Rq 
R\z=o , and f R0 = fii(Ro)- We further infer 



C l - 2 o 2 - 

— m — 2k pa 

C2 



(13) 



from requiring equivalence with ACDM when \f R o\ ~> 0. 
From Eq. ((7]) , we get 



v(4>) = 



Rf R (l + l/n) + 2 K 2 p A 
2k 2 (1 + f R ) 2 



(14) 



V6k(1 + JrY 



where R/R Q = (W/fl) 1/( " +1) - For f R < 1 and sub- 
tracting the background, Eq. (O becomes 



V 2 $fn = 3 [«*(/*) - « 2 fcm] 



(16) 



where Sf R = f R (R)-f R (R), 5R = R-R, Sp m = p m -p m . 
We assume a spatially homogeneous and isotropic cos- 



mological background metric and consider perturbations 
of the Friedmann-Lemaitre-Robertson- Walker line ele- 
ment in longitudinal gauge, i.e., ^ = Sgoo/(2goo) and 
$ = 8ga/(2gii). Combining the time-time and time- 
space component of the perturbed modified Einstein 
equations, one obtains 



J 



3f' R * + 2V 2 [{l + f R )V]-3Hf' R y = K 2 (Sp + 3Hv) 



QH 2 



3/r 



(l + /«) 2 



SfR 



S(f R R-f) , 3f R 



r 



i + /i 



(17) 



We assume matter dominance, i.e., dp = Sp m and v = v m , 
the matter peculiar velocity potential defined by d^v = 
—dun, where u^ is the unit four- velocity. Here, dots de- 
note derivatives with respect to cosmic time. For /jj<1 
and in the quasistatic limit, V 2 3> H 2 , neglecting time- 
derivatives of the perturbations and f R , this yields the 
modified Poisson equation 



V 2 * = ^5p m -l5R 



R) 



(18) 

where we have used Sf ps f R 5R and Eq. (|H)j) to replace 
V 2 <%. 

Note that if the background field \fm\ is large com- 
pared to typical gravitational potentials (~ 10~ 5 ), we 
may linearize the field equations, Eqs. (|i"6|) and (|T5|) . via 
the approximation 

5 /.R R=R 

where m is the mass of the f R field at the background. In 
Fourier space, the solution to Eqs. (|T6"|) and (|T5|) within 
the linearized approximation is 



<5i?: 



5/r = 3m 2 Sf F 



(19) 
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1 



a 2 <5p m (k), 



(20) 

for a comoving wavenumber k = |k|. For scales k 3> ma, 
this leads to an enhancement of gravitational forces by 
a factor of 4/3. Computations using Eq. (|2U)) are re- 
ferred to as the no-chameleon or linearized f(R) case [161 ] . 
whereas solutions to Eqs. (|16|) and (|T8|) are referred to as 
chameleon f(R) gravity. In the following section, we will 
study solutions for the scalar field Sf R and the Newtonian 
potential "J within a virialized cluster in both scenarios. 



III. CLUSTER QUANTITIES 

Effects from f(R) modifications of gravity on halo 
properties were studied in, e.g., |19l - (23j . The enhanced 



abundance of clusters caused by the modification was 
used in [LJ, |l4| to place an upper bound on the scalaron 
background value of |/ro| ;$ 1CP 4 . Ref. [12| used cluster- 
galaxy lensing measurements of the excess surface mass 
density to constrain \f R o\ based on the f(R) enhance- 
ment of the cluster density profile around the virial ra- 
dius (see |20|), finding an upper bound of \fm\ < 10~ 3 . 
These analyses have been carried out in the linearized 
regime of f(R) gravity, i.e., where the approximation 
Eq. (|19j) is valid. However, with future measurements, 
constraints will penetrate into the chameleon regime and 
it becomes important to incorporate the effects of the 
chameleon mechanism on the obscrvablcs. This shall be 
the concern of this section. Based on the NFW profile, we 
derive here analytic formulae, e.g., for the gravitational 
potential and velocity dispersion as observed within the 
virialized cluster. These relations can subsequently be 
used to predict observables in f(R) gravity without hav- 
ing to rely directly on simulations and constrain \f R o\ in 
the chameleon regime without the necessity of running a 
large number of simulations. 



A. Density 

Navarro, Frenk, and White (NFW) found that the 
simple relation 



Sp m {r) 



P Pc 



X-(l + Z. 



(21) 



provides good fits to the cluster density profiles measured 
in Newtonian CDM simulations. Here, r s denotes the 
characteristic scale and p s = j3 p c is the characteristic 
density with the critical background density p c . As we 
will show in ^IVl this simple function provides compara- 
bly good fits to f{R) gravity simulations in both the lin- 
earized and chameleon scenarios for r G (ro,r v ; r ), where 
r v ; r is the virial radius of the cluster and tq is conser- 
vatively set by the requirement that N > 800 particles 



in the simulations fall within that radius (see £HVI and 
Appendix [B]) . In Appendix [SJ we argue that the appli- 
cability of the NFW profile to modified gravity can be 
motivated by consideration of the self-similar infall of a 
collisional gas under modified forces, preserving consis- 
tency with the Jeans equation and the virial theorem. 

In the following, we shall assume this profile to be 
exact on the scales of interest, r £ (ro,r v ir) ; and de- 
rive from that (semi-) analytic relations for the mass, the 
scalar field, the gravitational potential, and the velocity 
dispersion. 



B. Mass 

We integrate Eq. (pH]) from the origin of the cluster to 
the radius r to obtain the mass of the cluster enclosed by 
r' < r. However, the NFW profile might not apply for 
r' < ?*o. For fair comparison with simulation results, we 
add a correction for the inner part of the cluster. The 
mass can then be obtained by the integration 

M(r) = 4tt / dr'r' 2 [S Pm (r') + A Pm (r')e(r a -r')} 
Jo 



as well as f(R) gravity. 



47r/3p c r s 3 



In 1 



+ M C , (22) 



where Ap m (r) = Sp mtSim (r)-Sp m ^ FW (r) and M c defines 
a mass calibration at tq. Eq. (|22[) applies to Newtonian 



Scalar field 



We first derive the scalar field 5 fa in the linearized 
case and based on this result, we then construct an an- 
alytic approximation for the solution in the chameleon 
case by requiring that Sfn 7> |/ro|- For comparison, 
we also study the approximation of the chameleon scalar 
field of Pourhasan et al. [25| and a numerical solution 
to the scalar field equation, Eq. (fT6|) , assuming spherical 
symmetry and the applicability of the NFW halo den- 
sity profile. For clarity, we shall denote our solutions for 
the scalar field as <5/^ n and Sf^ am for the linearized and 
chameleon case, respectively. 



1. Linearized field 

In order to find <5/]j n , we solve Eq. (|16[) with the ap- 
proximation Eq. (|19|) and the assumption that Sp m is de- 
scribed by a NFW profile. Furthermore assuming spher- 
ical symmetry, we obtain the differential equation 



dl + -d r )5fi?-m 2 5f£ + 



P* 



3 r(r + r s ) 2 



0, (23) 



which has the solution 



J 



SfA n (r) 



^p s r 3 s {r [0, m(r + r„)] e 2m ^ +r ^ + T [0, -m(r + r„)] } + (Ci + ^ 



G 



2 2m r \ m r 



D -m(r+r a ) 



r 



(24) 



Here, C± and C2 are integration constants and 

/*OQ 

r(s,r) = / dtt s 



LS-l e -t 



(25) 



denotes the upper incomplete gamma function. Ci is 
the amplitude of a growing mode and since we want to 
restore general relativity (GR) at r — > oo, we set C-2 = 0. 
We further require K 2 5p m /3 to dominate over rn 2 5f^ 1 
towards the origin of the halo and hence, 



(26) 



limr<5/]j n ->•(). 



This is also apparent from the resemblance of Eq. 
when m 2 5fft n is small to the standard Poisson equa- 
tion and its solution for the Newtonian potential ^gr 
(see EIIIIDI) . The integration constant then becomes 

2 3 

Cx = -!-PfjL e -™r s [ e 2m ^r(0,mr s ) + r(0,-mr s )] 



G 



(27) 



and hence, we arrive at our solution for the linearized 
scalar field 



SfA n (r) 



2 3 

ft p s r s 
6 



{r[0,rn(r + r s )]e 2m(r+ra) 



+r[0, —m(r + r s )] — T(0, —mr s ) 



2m r g 



r(0, mr s )} 



-?7i(r+r s ) 



(28) 



Note that at scales where p B 
obtain the limits 

8f%> - 



p m and p m > p m , we 



3m 2 

„2 



8pn 



V 2 5f™ = ~5p m , 
respectively Eq. (J3UJ) implies that 

r ,iin _ ^ 2 p s r s 3 ln(l + r/r s ) K\ 
6JR ~~3 r + V 



(29) 
(30) 

(31) 
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FIG. 1: Properties of the scalar field 8fn for a cluster of M v i r = 1-36 x I0 14 M^/h and a background field amplitude of 
|/flo| — 10~ 5 (long-dashed line) with n — 1. Le/t: The analytically derived scalar field 5 fit (solid line) and its behavior at the 
limit of large and small scales, respectively. The transition of the linearized field into a chameleon field occurs instantaneously 
where Sfn = \fno\- Middle: Comparison between the instantaneous chameleon transition (solid line) and the (^-transition for 
8fn (dashed line) with analog semi-analytical derivation to [25[ • Matching 5fn = 5/.R,sim at r vlr brings the approximations into 
good agreement with the simulated 5fn(r) (dots). The thick long-dashed lines show numerical solutions to Eq. (|16[) assuming 
a NFW profile with integration constraints set by additionally requiring that 5f' R (r v i r ) corresponds to the analytic and semi- 
analytic result, respectively. Right: Enhanced force by the linearized (dotted line), the instantaneous chameleon (solid line), 
and the C 1 (dashed line) chameleon scalar field <5/r (matched at TVir), respectively. The lines of the linearized and instantaneous 
chameleon field overlap beyond the chameleon region. 



where the integration constants Ki = due to Eq. ([25)1 
and i\2 = — m K 2 p s rg exp(?rir s )r(0, r7ir s )/3 to match 
Eq. ([25)1 at the origin. Hence, for p m ^S> p m , 



*m 



tin 



3 



*GR- 



« 2 Psrf 



TOe m,a r(0,TOr s ), (32) 



where ^gr is taken for an isolated halo assuming an 
NFW profile on all scales (see 3111 Dp . 

We illustrate our solution for the scalar field and its be- 
havior at the limit of large and small scales, respectively, 
in Fig. HJ 



2. Chameleon field 

In order to describe the chameleon mechanism in f(R) 
gravity, let us revisit Eq. (fTU| . which for |/jj| <C 1 and 
a = -k/V6 becomes 



= -— Pm + V'(t). 



(33) 



In high density regions, where k p m ^> — \/6V 2 </>, the 
scalar field becomes 



jj K = /ho 



Ro 



K 2 (p m + 4p A ) 



n+1 



Jr- 



(34) 



Since p m ^> p c , for n > — 1, we get f R ~ K<p — 0. Hence, 
modifications of gravity are suppressed. 

We consider three different approaches for describing 
the transition of <5/]^ n to Sf R ham and compare them with 
each other. The first approach is the assumption of a sim- 
plified, but analytically describable, instantaneous tran- 
sition to Sfu = — fm, whereas the second is a semi- 
analytical match of the two regions and the third is a 



numerical solution to Eq. (|16[) . Thus, in the first case, in 
order to implement the chameleon mechanism in our fit, 
we may simply require 



/■el 
JR 



chain 



(/iM 



(35) 



or equivalents, 6 ft™* = min (5f™, -f R0 ) (see Fig. [J). 
As we will show in 3IV1 this yields a good approximation 
to the simulated chameleon f R field. The difference being 
a more gradual decrease in d r f R from — 29 i - , J/gr/3 to 
instead of an instantaneous transition (cf. Fig. Q]). The 
chameleon transition is, however, very efficient, which 
allows the applicability of this approach. 

For comparison, we follow the semi-analytic approach 
of Pourhasan et al. [25[ for describing the chameleon field 
for an inverse power-law potential of a scalar-field. Their 
procedure corresponds to matching the f(R) chameleon 
interior solution, Eq. (|34j) . applying to r G (r_, r+) to the 
chameleon exterior solution, Eq. (|3ip . for r > r c at the 
transition scale r c . This defines the integration constants 
in Eq. ([3T|) . K^ = follows from the requirement that 
Sfn — > when r — > oo. K\ and r c are determined by 
requiring that 



Sf R (r) = (Sf 



R,r<r, 



u^>, 



)(r)eC 1 (U) (36) 



with r c e U C Kq", i.e., the matched scalar field and its 
derivative are continuous at the transition. We compute 
r c numerically and show the according solution for dfu 
in Fig. [TJ Note that Eq. (|36|) assumes that the interior 
solution also applies to the shell r <E (r + ,r c ), where r + < 
r c . Following [25j , the boundaries of the interior region, 
i.e., the regime of applicability of Eq. (f34]> . r £ (r_,r + ), 
can be obtained from the roots of — 3V 2 /fl/K 2 p m = 10~ 2 . 
If for r + < r c , we have (r c — r+) <C r c , the interior 
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FIG. 2: Same as middle panel of Fig. [T] but with cosmological 
background density as boundary condition instead of match- 
ing to <5/fl,aim(»"vir) for 5/r. For the radial derivatives, we 
further require that 5f' R (r\i r ) corresponds to the analytic and 
semi-analytic result, respectively. 



solution may be extended into the shell, which in this 
case is sufficiently thin. Strictly speaking, the condition 
r + > r c is not satisfied for the scalar field shown in Fig. [1] 
The Sfu computed with this procedure, however, still 
yields a good description to the simulated scalar field 
(see Fig. [TJ. 

As a third case, we numerically solve the differential 
equation for S/r, Eq. (|16[) . assuming spherical symmetry 
and that 5p m is given by a NFW profile. For stability 
reasons we use the substitution /# = — e u ^ (see |17|) in 
our computations. We compare the numerical solutions 
for dfn obtained in this way to the chameleon scalar 
fields obtained through the analytic and semi-analytic 
approach, Eq. ([55]) and Eq. (|5tj|). respectively, in Fig.Q] 

Note that due to the limited applicability of the NFW 
density profile for the description of Sp m beyond r € 
(jo, r v j r ) and the unknown environment, when comparing 
to simulations, we correct the analytic and semi-analytic 
solutions of fa to match them to /_R,sim at r v j r , or equiv- 
alcntly, we require this constraint when numerically solv- 
ing Eq. (|16j) . In specific, this means adding a constant 
deviation from the background density in Eq. (|23[) and 
dropping the condition Ki = in the instantaneous and 
C 1 chameleon solution, respectively. This brings the com- 
puted scalar fields into good agreement with the simu- 
lated Sfn over the radial range r £ (rn, r v ; r ) as is demon- 
strated in Fig. [TJ For the comparison to the numeri- 
cal solution of Eq. (|16|) . we further constrain the inte- 
gration by requiring that 5f' R (r v i r ) corresponds to the 
derivative of the solution for the C 1 and the instantaneous 
chameleon field, i.e., the semi-analytic and analytic <5/ij, 
respectively. 

It is important to note that matching 5fn(r v i T ) to sim- 
ulations is essential for recovering the radial dependence 
of the scalar field from simulations. As demonstrated 
in Fig. [21 if we use the cosmological background as the 
boundary condition instead, i.e., Ki = in Eq. (|31[) , we 
arc not able to reproduce the scalar field <5/_R jS im- 



D. Gravitational potential 

For a spherically symmetric mass distribution, the 
gravitational potential at r £ (^Oi^vir) in Newtonian 
gravity is obtained by the sum of the interior and ex- 
terior spherical mass shells, i.e., 



*GR 



dM(r') - 



dM(r') 



(37) 



dr> r' 2 [p m (r>) + Ap m (r')6(r - r')} 



dr' r' [p m (r') + Ap m (r')6(r' - r vir )] }>(38) 
K 2 p s r^ln(l + r/r s ) k 2 M c 



8tt r 



*c + *ext. (39) 



Here, 'Fcxt indicates an external gravitational field in case 
the halo is not isolated and ^ c accounts for deviations 
from the NFW density profile at r > r v - u ., e.g., the two- 
halo contribution, i.e., \& c does not vanish even if the halo 
is isolated. 

Combining Eqs. (JTSJ) and (JTHJl yields 



V 2 * = U Pm l -V 2 5f R 



= V* 



*GR - -6f R 



(40) 
(41) 



By partial integration and with the analog boundary con- 
ditions as in the integration of 'FgRj i-e., 

]imrd r Sf R = 0, (42) 

lim d r 5f R - 0, (43) 



we obtain the modified gravitational potential 

* = *gr - 2*/r. (44) 

For comparison to simulations, the combination , F C + 
"Foxt, i.e., the halo density correction to the NFW profile 
beyond r v ; r and the external gravitational field 'Fext from 
the environment, is calibrated to 'Fsim(rVir)- 



E. Velocity dispersion 

From contemplations on the self-similar secondary in- 
fall and the shocked accretion of a collisional gas onto 
the center of an initially spherical uniform overdensity 
in an otherwise unifor mly expanding Einstein-de Sitter 
universe, Bcrtschinger [26| determined power-law behav- 
iors for the fluid variables, which produce the phase-space 
density profile 



P 



5/2 



,3/2 



-15/8 



(45) 



where p denotes the pressure of the gas. We refer the 
reader to Appendix fX] for details about this derivation 
and its applicability to modified gravity. With a gas of 
pressure p = p m a 2 , where a is the velocity dispersion, 
Eq. (|45|) yields the pseudo phase-space density (PPSD) 
profile 



a(r) s 



lA 



±£± (H) 

4 of V r ) 



15/8 



(46) 



with a s = cr(r s ). The simple relation Eq. (J46J) was found 
to give good fits to Newtonian CDM simulations |27j . 

As pointed out in Appendix [AJ the r-dependence of the 
PPSD profile does not change in a simplified approach of 
f(R) gravity, i.e., where a simple force amplification of 
4/3 is assumed. We therefore use the velocity dispersion 



I 4p m 

Ps 



2/3 



5/4 



(47) 



to compare to simulations, where a s is taken to be an 
additional degree of freedom that we fit to simulations 
and /o m is assumed to be correctly described by the NFW 
density profile for r G (ro, r v - 1T )- 

Note that the radial effect on a 1 from a transition 
in the modified forces is blurred out over a wide range 
of scales. As shown in mVl the radial dependence of 
Eq. (|47|) fits the chameleon f(R), linearized f(R), and 
ACDM simulations equally well. For the estimation of 
cr s , however, in order to encompass the chameleon as a 
function of halo mass, we can replace the constant force 
enhancement with a weighted average over the modifica- 
tion shown in the right panel of Fig. Q] (cf. |21|). 



A. iV-body simulations 

The simulations used in this work were carried 
out for the Newtonian (GR), linearized (N), and full 
chameleon (F) scenarios for each field strength \Jrq\ = 
fO- 6 ,10- 5 ,lCT 4 with n = 1 [l7|. Each set of simu- 
lations consists of 10 realizations with each box size, 
ibox = 64ft,- 1 Mpc, 128/1- 1 Mpc, 256/i _1 Mpc, and a 
total particle number of N p = 256 3 placed on 128 3 do- 
main grids. Thereafter, the different set of simulations 
arc denoted by GR-[L box ], N-[-log 10 |/jjo|H-kbox], and 
F-[— log 10 1/ijolH-kbox]- During the simulation, the do- 
main grids arc refined progressively where the local densi- 
ties are sufficiently large to reach a predefined threshold. 
In this way, the grid structure efficiently follows the den- 
sity distribution so that the high density regions can be 
well resolved. The cosmological parameters are fixed to 
values following the WMAP 3-year results, £7a = 0.76, 
f2 m = 1 — £1a ; h = 0.73, n s = 0.958, and the initial 
power in curvature fluctuations A s = (4.89 x 10~ 5 ) 2 at 
A: = 0.05 Mpc" 1 . 

Halos within the simulation and their associated 
masses are identified via a spherical overdensity (SO) al- 
gorithm (cf . [28[ ) . The particles are placed on the grid by 
a cloud-in-cell interpolation and counted within a grow- 
ing sphere around the center of mass until the required 
overdensity is reached. The mass of the halo is then 
defined by the sum of the particle masses contained in 
the sphere. This process is started at the highest over- 
density grid point and hierarchically continued to lower 
overdensity grid points until all halos are identified. Note 
that the virial overdensity obtained for ACDM is used to 
identify halos even in f(R) gravity in order to make a 
fair comparison between the different models (see Ap- 
pendix [Bj . 



IV. COMPARISON TO SIMULATIONS 

Based on the NFW fit for the halo density profile, 
Eq. (|21l) , we have constructed in i]IIII analytic fits to the 
halo mass enclosed at radius r, Eq. (|2"2"j) , the scalar field 
Sfn, Eqs. ([28]) and ([35]) . and the modified gravitational 
potential given by Eqs. ([55)1 and (|4"4")l. We have further 
assumed that the velocity dispersion of dark matter parti- 
cles is correctly described by the power-law PPSD profile 
predicted by the self-similar collapse of a collisional gas. 

In this section, we shall test these relations against 
dark matter TV-body simulations of Newtonian, linearized 
f(R), and chameleon f(R) gravity. Thereby, we assume 
the overdensity [3 (or equivalently p s ) and the character- 
istic scale r s in the NFW profile, as well as the velocity 
dispersion at r s , i.e., c s , to be free fitting parameters. 
We then compare the quality of these fits between the 
different simulation outputs and analyze the ability of 
scaling relations based on the spherical collapse to give 
predictions for p s , r s , and a s . 



B. Performance of fits 

We test the predictions made in EjIIII on the z = 
simulation output described in £ )IVAI We calculate the 
reduced \ 2 of the relative deviation between the predic- 
tion for the quantity q from our analytic fits and the 
simulation output, i.e., 



Ag,rcd 



1 N / \' 



(48) 



for each halo of the simulation independently, where N 
is the number of bins in r € (ro, r v ir) and v = N — n — 1, 
n being the number of fitting parameters used in the fit. 
We then calculate the mean of Eq. (f4"5)) over all halos. 
Our results are summarized in Table [U The mass range 
chosen for the selection of the clusters (see Table H| picks 
about 40-50 halos out of the 50-60 most massive halos 
of each simulation. The average is then taken over the 
10 realizations of each simulation configuration, leading 
to an average over about 500 halos for the results shown 



Simulation 






Mass 


range 


(M) , 


J \X-6p m ,red) ' 


J \Xm, red/ ' 


/far-) ' 


\J \ X Sf R ,rcd) • 


J \Xy, Ic d/ 


GR-64 


2 


X 


10 13 


-2 


X 


10 14 


4.25 


X 


10 13 


10 


17 


10 






GR-128 


1 


X 


10 14 


-5 


X 


10 14 


1.74 


X 


10 14 


8 


14 


9 






GR-256 


3 


X 


10 14 


-7 


X 


10 14 


4.14 


/ 


10 14 


5 


9 


9 






N-64-4 


2 


X 


10 13 


-2 


X 


10 14 


4.44 


X 


10 13 


10 


18 


10 


1 


6 


N- 128-4 


1 


X 


10 14 


-5 


X 


10 14 


1.84 


X 


10 14 


8 


15 


9 


1 


5 


N-256-4 


3 


X 


10 14 


-7 


X 


10 14 


4.21 


X 


10 14 


5 


10 


8 


1 


3 


N-64-5 


2 


X 


10 13 


-2 


X 


10 14 


4.35 


X 


10 13 


10 


17 


9 


4 


6 


N- 128-5 


1 


X 


10 14 


-5 


X 


10 14 


1.82 


X 


10 14 


8 


15 


9 


4 


5 


N-256-5 


3 


X 


10 14 


-7 


X 


10 14 


4.19 


X 


10 14 


5 


10 


9 


4 


4 


N-64-6 


2 


X 


10 13 


-2 


X 


10 14 


4.45 


X 


10 13 


10 


17 


10 


7 


6 
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1 


X 
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-5 


X 


10 14 


1.76 


X 
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8 


14 


10 


7 


6 


N-256-6 


3 


X 
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-7 


X 
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4.14 


/ 
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5 


9 


10 


5 


3 
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2 
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-2 


X 
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10 


6 


6 
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-5 
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10 14 


10 


14 


8 


5 


5 


F-256-4 


3 
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-7 


X 
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4.18 
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5 


10 


8 


4 


3 


F-64-5 


2 


X 
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-2 


X 


10 14 


4.33 


X 
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10 


18 


9 


4 


6 


F-128-5 


1 


X 


10 14 


-5 


X 


10 14 


1.71 


X 


10 14 


7 


15 


7 


4 


4 


F-256-5 


3 


X 


10 14 


-7 


X 


10 14 


4.15 
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10 14 


5 


10 


7 


2 


3 


F-64-6 


2 


X 


10 13 


-2 


X 


10 14 


4.38 


X 


10 13 


10 


16 


8 


1 


5 


F-128-6 


1 


X 


10 14 


-5 


X 


10 14 


1.73 


X 


10 14 


8 


14 


8 


0.1 


4 


F-256-6 


3 


X 


10 14 


-7 


X 


10 14 


4.12 


X 


10 14 


5 


9 


9 


0.1 


4 



TABLE I: Comparison of goodness of fit of the analytic predictions of £11111 based on the NFW halo density profile and the 
power-law PPSD between the GR, linearized f(R), and full chameleon f(R) simulations. Cluster masses are given in M^jh. 
\Jx\cd i s computed from the %-deviation from simulations for the halo profile 5p m , the mass M, the velocity dispersion a, the 
scalar field 5 fa, and the gravitational potential $ inr £ (ro,r v ir)- 



in each row of Table Q] We chose this simple approach 
of quantifying the goodness of fit only for the qualitative 
comparison between its performance in the concordance 
model and in the linearized and chameleon f(R) gravity 
models. 

In particular, we find that both the NFW profile and 
the power-law PPSD profile with standard radial depen- 
dence, r -15 / 8 , yield equally good descriptions to the ha- 
los produced in the f(R) models as they do in the con- 
cordance model. Note that for the scalar field Sfu, we 
only take into account the instantaneous transition to the 
chameleon region (see SflTTC]) with the mass calibration 
at ro , Eq. (f2"2"j) , in correspondence to the mass correction 
in the gravitational potential Eq. ([55)1 . 

For illustration of our fits on the simulation data, we 
choose a set of simulations that highlights the effect 
of the chameleon mechanism, i.e., where the transition 
from the large-field to the small-field limit takes place 
within the scales of interest. Thus, we seek a combi- 
nation of medium field strength |/ro| and medium halo 
masses. We therefore illustrate the halo quantities pro- 
duced for \Jro\ = 10~ 5 and Lb ox = 128 Mpc/h, cor- 
responding to the set of simulation outputs denoted by 
F-128-5 (see gV"A"|) . Fig.|3]shows the stacked fits to the 
F-128-5 simulation output, which is also stacked, for the 
halo density profile, the cluster mass, the velocity disper- 



sion, the PPSD profile, the scalar field fa, and the grav- 
itational potential. The normalizations and narrow mass 
range, (1.65 — 1.70) x 10 14 M & /h, are chosen such that 
standard deviations of the simulation output are small, in 
particular, at scales where the chameleon mechanism is 
active. The narrow mass range is applied to the 10 real- 
izations of the F-128-5 configuration, resulting in taking 
the average over 16 halos. 

In Fig. 2] we show the overdensity Sp m /p m at r s , the 
characteristic scale r s , the velocity dispersion squared at 



07, and the concentration parameter c v 



r/r s 



as a function of the virial mass M v j r . The three types 
of data points correspond to the Newtonian, linearized 
f(R), and full chameleon f(R) model, respectively. The 
three bins of each type indicate the stacked best-fit val- 
ues to each halo produced in the simulations of the three 
different box sizes, where the least massive ones corre- 
spond to Lbox = 64 Mpc/h and the most massive ones 
to £box = 256 Mpc/h. We stack the about 10 most mas- 
sive halos in the selection for Table Q] of each realization 
of each simulation configuration, resulting in the average 
over about 100 halos for each data point. Fig. 0] is a 
good demonstration of the chameleon mechanism in sev- 
eral realizations. For the f(R) models, we clearly observe 
a shift of the average mass of each data point towards 
higher masses with respect to GR. This corresponds to 
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FIG. 3: Comparison of radial dependencies of f(R) halo quantities predicted by simulations and the fits constructed in £11111 
for | /ho | = 10~° (n = 1) and Lbox = 128 Mpc/h (F-128-5). The simulated halos and their individual fits are stacked for 
M = (1.65 — 1.70) Mq//i. The error bars show the standard deviation in the simulation output. Top left: cluster density 
profile Spra/pm- Top middle: halo mass M. Top right: velocity dispersion a. Bottom left: pseudo-phase-space distribution 
p/a 3 normalized at r B . Bottom middle: fn scalar field for the stacked instantaneous (solid) and C 1 (dot-dashed) transition to 
the chameleon solution. Bottom right: gravitational potential \& for the 6/r solutions in the bottom middle panel along with 
the limiting assumptions of Newtonian/GR (dotted) and linearized f(R) (dashed) forces. 



the enhanced abundance of massive halos in f(R) grav- 
ity. Whereas for high values of |/ro| the full chameleon 
simulations approach the linearized simulations, they re- 
produce the Newtonian simulations at low values of \fm\. 
We further observe that in the full chameleon simula- 
tions, the displacement in the mean mass with respect 
to GR is strongest at high masses for large \fm\ and at 
low masses for small \Jro\, respectively. This coincides 
with the expectations of the f(R) halo abundance, i.e., 
where chameleon simulations recover GR at high masses 
for low values of \fm\ but differ at low masses; and for 
large values of |/i?o|, the strongest difference to GR is ob- 
served at high halo masses (see [13, [2CJ ) • A further real- 
ization of the chameleon effect can be observed in the ve- 
locity dispersion. Here, the predictions of the chameleon 
simulations correspond to the linearized simulations for 
| /ho | = 10~ 4 that are enhanced by a factor of ~ 4/3 
over the GR predictions. For \/rq\ < 10 -5 , however, 
the chameleon simulations recover the GR velocity dis- 
persion as the f(R) modification is suppressed and the 
gravitational force returns to being Newtonian. 



C. Prediction of fitting parameters 

We predict the virial halo concentration c V h- along with 
the two fitting parameters p s and r s in f(R) gravity and 
GR using scaling relations based on spherical collapse 
calculations (see Appendix [B]) following [2fJ . In order to 
determine the third fitting parameter, the velocity dis- 
persion at the characteristic scale r s , o~ s , we require that 
the Jeans equation must be satisfied at r s given the as- 
sumptions of a NFW halo density profile and the stan- 
dard radial dependence of the PPSD profile, along with 
a fit for the velocity dispersion anisotropy relation based 
on concordance model simulations (sec Appendix [5]) . 

We find that the scaling relations obtained in this way 
yield qualitatively good reproductions of the best-fit val- 
ues of p s , r s , c s , and c v ; r . Our results are shown in Fig. 2] 



V. CONCLUSIONS 

Modifications of gravity have extensively been tested 
on solar-system scales (see, e.g., |20|) and to a lesser de- 
gree at large cosmological scales using specific alternative 
theories of gravity (e.g., j!3l Il4l . l3Q -- [32j ) . as well as generic 
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FIG. 4: Stacked halo fitting parameters along with the concentration c v i r for the three different box sizes and predictions for 
them from scaling relations (see Appendix [Bj . In f(R) gravity, halo abundances and therefore mean masses are enhanced. 
Clearly identifiable, in the linearized f(R) regime, velocity dispersions are enhanced by a factor of 4/3 over the Newtonian 
results, i.e., in addition to smaller effects of different c v i r . The chameleon effect returns the cluster abundance at high halo 
masses and the velocity dispersion to the Newtonian/GR predictions. The predictions from scaling relations based on the 
spherical collapse and the Jeans function are shown for the Newtonian (GR - solid line) and linearized f(R) (N - dashed line) 
case, respectively. Rows from top: The matter overdensity at r s (given by /3 or p a ), the characteristic scale r a , and the velocity 
dispersion squared at r B , of, as well as the concentration c v i r - Columns from left: \fso\ — 10~ 4 , 10~ 5 , 10 -6 (n = 1). 
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modifications to GR while adopting a ACDM background 
(e.g., |33l - [38| ) or simultaneously allowing a dynamic ef- 
fective dark energy equation of state [33, E(| . However, 
gravity may also be tested by the structure observed at 
intermediate scales [IJ, [l2], |4J| . In this regime, nonlinear 
gravitational interactions gain in importance and need 
to be modeled correctly to obtain reliable predictions for 
both GR and its competitors, which in turn can be com- 
pared with observations to infer constraints on modified 
gravity theories. Besides the usual difficulties of model- 
ing the nonlinear structure known to studies of the stan- 
dard Newtonian gravity, modifications of gravity may be 
complicated by additional nonlinear mechanisms such as 
the chameleon effect, which suppresses modifications of 
gravity in high density regions. In order to obtain reli- 
able constraints from observations, these effects need to 
be consistently incorporated. 

In this paper, we concentrate on f{R) gravity and aim 
at describing the scalar field, the gravitational poten- 
tial, and the velocity dispersion within virialized clus- 
ters produced in f(R) JV-body dark matter simulations 
of the Hu-Sawicki model. We derive and test analytic, 
semi-analytic, and numerical relations for these quanti- 
ties that can be used to compare theory with observations 
at the virialized scales of clusters. We assume the stan- 
dard NFW halo density profile and the PPSD profile with 
usual power law, which we find to provide comparably 
good fits to the f(R) scenario as they do for the concor- 
dance model. We argue that this is not unexpected from 
the consideration of modified forces in the secondary in- 
fall of a collisional gas describing the TV-body dark matter 
system. The fits to the simulation output are based on 
three degrees of freedom, the characteristic density p s , 
the characteristic scale r a , and the velocity dispersion at 
r s , a s . We find that scaling relations based on the grav- 
itational collapse and the requirement of the validity of 
the Jeans equation yield good qualitative predictions for 
these quantities when accounting for the modified forces 
at work. 

The extension of our results to scales beyond the virial 
radius exhibits additional challenges, such as the correct 
modeling of the two- halo contribution (cf. [12|,[20[), which 
we shall leave for future work. 



Appendix A: Density profiles in modified gravity 
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In i jIVBl wc have seen that the NFW halo density 
profile for r £ (r , r v j r ) provides as good fits to the virial- 
ized clusters in f(R) gravity as it does in the Newtonian 
scenario. In the following, we shall give qualitative argu- 
ments for why this may be expected from a theoretical 
point of view. In order to study the gravitational col- 
lapse in f(R) gravity, we parametrize the enhancement 
of gravitational forces by the factor (1 — F), where for 
simplicity, the quantity F is defined by 



F = 



-1/3, /(i?) gravity, 
0, Newton/GR, 



(Al) 



which holds in the linear regime and when modifications 
are absent or suppressed, respectively (see Fig. [IJ. 



1. Self-similar infall of a collisional gas 

We follow Bertschinger |26[ for the self-similar infall 
and the shocked accretion of a collisional gas onto the 
center of an initially spherical uniform overdensity in 
an otherwise uniformly expanding Einstein-de Sitter uni- 
verse. The equations governing the motion of the fluid 
are, nondimensionalized (see |26[), 



(V - £a)£>' +DV+ =^— - 2D 
9 A 

Q 1 p/ 



0, 
2M , 



(V--X) 



T 



7- 



iy_ 
15 
w 



20 
3A 2 D 



27, 



(A2) 



i.e., the continuity, Euler, adiabatic, and mass equations. 
Here, V, D, P, and M is the nondimensionalized velocity, 
density, pressure, and mass, respectively, and 7 indicates 
the ratio of specific heats, which is taken to be 7 = 5/3, 
i.e., the ratio for a monatomic ideal gas. Primes denote 
derivatives with respect to A = r/r ta , where r ta is the 
turn-around radius. Note that we have introduced here 
the modification of the gravitational force F. The nondi- 
mensionalized quantities relate to the fluid variables via 
(see dj) 



v(r,t) 



rta 



^(A), 



p(r,t) = PEdsD(X), 



p(r,t) 



PEdS 
47T 



(t) f «' 



m(r,t) -- yp EdS r t 3 a A/(A), 



(A3) 



where PEdS = 4k 2 t 2 /3 is the Einstein-de Sitter back- 
ground density and t is the cosmic time. 
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Factoring out the asymptotic behavior of Eqs. 
the origin, requiring the boundary conditions 

V = M = 0, A = 0, 



2) at 



(A4) 



it is easy to see that 
D = \- 9/4 D(\), P = A- 5/2 P(A), 



M = A 3/4 M(A) 
(A5) 
with finite D, P, and M at A = 0. Note that chang- 
ing F, changes the relation between P and D too, but 
dependencies on A remain the same, i.e., the nondimen- 
sionalized radial dependence of the density profile is not 
affected by the force modification. Therefore, for equal 
nondimensionalized mass, 



D = D. 



GR, 



P = P GK (1-F). 



We compare halos in the different gravitational models, 
however, by equating the corresponding virial (dimen- 
sional) masses with each other. But since we define the 
virial radius by the same virial overdensity A v ; r in all of 
the models (see Appendix [B]) , assuming equivalent back- 
ground, the virial radius r v ; r and therefore r ta are the 
same. For a gas of pressure p — per 2 , this therefore im- 
plies 



p{r) = pGR.(r), cr(r) 2 = a GR (r) 2 (l - F). 



(A7) 



Note, however, that Eq. (|A5|) implies p ~ r -9 / 4 , which 
does not provide a good fit to simulated dark matter ha- 
los, but we shall assume for now that the relations in 
Eq. (|A7|) hold even in cases where p(r) is not described 
by a simple power law. As pointed out in Jj A 21 and £|A31 
this assumption remains consistent with the Jeans equa- 
tion and virial theorem, respectively. Rather than the 
directly predicted radial dependence in Eq. (|A5I) . we are 
interested in the relation D 1 ^ 2 /P'^ 2 , which defines the 
PPSD p/a 3 . According to the relations in Eq. (|A5[) , 



P( r ) 
a(r) 3 



,-15/8 



(A8) 



which in turn is found to yield a good description for the 
results obtained from CDM simulations [27| . With the 
force modification F, we find that 



p(r) PgrW 



1 



a(r) 3 <J GR (r) 3 (l-F) 3 / 2 ' 



2. Jeans equation 



(A9) 



of energy-momentum and thus applies to all metric theo- 
ries of gravity and hence the f(R) model considered here. 
In spherical coordinates and with the force enhancement 
F, the Jeans equation can be written as 



Va 2 r = 
V = 



(1- 
d 
dr 



F) — y GR , 
dr 



dlnp 
dr 



2ft 

r 



m , (MO) 



where ay denotes the radial component of the velocity 
dispersion a, V defines a linear operator, and /3<r(r) de- 
notes the anisotropy in the velocity dispersion, i.e., 



1 



2a 2 



(All) 



(A6) Using Eq. (|A"9|) . we infer 



V 



{- 

\PGR 



2/3 



2 

"r.GR 



-^*GR = ©GR< GR) (A12) 



where we have divided the Jeans equation by (1 — F) and 
2^GR is t ne linear operator with p GR and /3 ct .gr- Hence, 
for p = pqr and a 2 = a 2 GR (l — F) , the Jeans equation 
is satisfied, implying /3 CT = ^ a ,GK- Note, however, that 
this solution is not the only one satisfying Eq. (IA10JI . 
Describing the alternative solutions is, however, beyond 
the scope of this paper. 



3. Virial theorem 

Multiplying the collisionless Boltzmann equation by 
velocity and position and integrating over both, we ob- 
tain for a system in steady state the virial theorem 



B^GR 



-2T ( 



GR, 



(A13) 



where Wqk and Tqr are the potential and kinetic ener- 
gies in Newtonian gravity, respectively, i.e., 

Wgk = - /"d 3 xp GR (x)x.V$ GR (x), (A14) 

r GR = iy'rf 3 xp GR (x) ( 7 2 JR (x). (A15) 

With p = pgr and a 2 = (1 — F )c GR , we have W = 
(1 - F)W GR and T = (1 - F)T GR , which consistently 
reproduces the virial theorem W = — IT . Following from 
the Boltzmann equation, this is as expected since using 
the Jeans equation in the integration of Eq. (|A14p or 
Eq. (|A15|) leads to the virial theorem, Eq. (|A13|) . 



The Jeans equation is derived from the collisionless 
Boltzmann equation, which describes the particle phase- 
space distribution as a function of position, momentum, 
and time. Thereby, the Boltzmann equation is multi- 
plied by and integrated over the velocity. The collision- 
less Boltzmann equation is the analog to the conservation 



Appendix B: Scaling relations from spherical 
collapse and the Jeans equation 

The three degrees of freedom used in the fits of £|III[ 
the amplitude p s (or /3) of the NFW halo density pro- 
file, Eq. ([2"TT) . the characteristic scale r B , and the velocity 
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dispersion at r s , er s , can be predicted by scaling rela- 
tions based on the spherical collapse (see, e.g., [20|]) and 
the Jeans equation. Here, we assume the usual collapse 
density S c ~ 1.673 and virial overdensity A v ; r w 390 
inferred from spherical collapse calculations with a stan- 
dard force, i.e., F = 0. See |2(j| for derivations of S c 
and A v ir in the modified spherical collapse with enhanced 
force F = -1/3. 

The variance is defined by 



(r) 



d 3 fc 

(2^)3 



W(kr) P L (fc), (Bl) 



where W is the Fourier transform of the real-space top- 
hat window function of radius R, i.e., 



W(kR) = 3 



sin(fc_R) cos(fc-R) 



L ( kR Y 



{kRf 



(B2) 



and Ph(k) is the linear matter power spectrum of the 
model. We determine F\,(k) by using the Eiscnstcin- 
Hu transfer function |42|, 14311 and the parametrized post- 
Friedmannian framework J44j for the designer f(R) grav- 
ity model to determine the approximate growth needed 
to obtain the power spectrum (cf. [14J])- The radius r in 
Eq. (|B1[) is defined by the cluster mass M enclosed by it 
through 



■(M) 



3M 

4irp lr 



1/3 



(B3) 



which defines er var (A-f). 

Given the virial overdensity A v ir and a virial mass 
M v i r , the virial radius r vn . of the dark matter cluster 
is determined by M v [ r = Att p m A vir r^ ij: /3. We further 
assume a concentration c vlT = r V i r /r s given by |45j 



r (M vir ,z = 0) = 9 



M* 



0.13 



(B4) 



requiring cr va r(-Af*) = 5 C . This describes the character- 
istic scale r s as a function of the virial mass M v { r . The 
amplitude p s of the NFW profile, Eq. (|2"Tj) . is then defined 
by the integration to M v - lr . 

The fits provided by this scaling relation are shown in 
Fig.|H Note that the relation Eq. (|B4[) was obtained from 
fitting to halos of mass (10 11 — 10 14 ) Mq/H produced in 
Newtonian CDM simulations |45| and the applicability 
to the halos in this study is therefore limited. Never- 
theless, this scaling relation is found to qualitatively re- 
produce the fitting parameters of the two extreme cases, 



i.e., ACDM and linearized f(R) gravity. Note that due to 
differences in the linear matter power spectrum between 
GR and the f(R) model, the characteristic scale r s and 
therefore p s for a given M V - 1T changes for f(R) modifica- 
tions with respect to the Newtonian results. Hence, the 
relation p(r) = pgr(t) (see §A1|) does not hold anymore. 
Note, however, also that assuming modified forces in the 
gravitational collapse changes S c and A v ; r and counter- 
acts the changes in p s , r s , and c v - n - (see [2(j]) to some ex- 
tent. For the illustration in Fig. 21 we assume standard 
forces in the derivation of S c and A v j r , i.e., F = 0. 

Finally, in order to obtain the velocity dispersion at 
r a , we revisit the Jeans equation, Eq. (|A10[) . for a clus- 
ter produced with Newtonian gravity. Assuming a NFW 
profile and the PPSD profile of Eq. (j46| . we obtain 



°"r,GR,s = ^PGR/s 



61n2-3 
25 - 24/3 CTiGR , s ' 



(B5) 



where /3 ct ,gr,s = /3o-,GR( r s)- We further assume that at 
r s , the velocity anisotropy relation is correctly described 
by0 



At,GR,i 



1 

40 



17 



23 din p 
6 dlnr 



37 
60' 



(B6) 



where dlnp/dln7' = — 2 at r = r s . Eq. (|B6[) was 
constructed as a fit to concordance model simulations 
in [46l | . The modified velocity dispersion is obtained 
from 0-2 s = (1 - F)a^ GRs with F = -1/3. We as- 
sume /3cr, s = /3 ffj GR.s and compare the predictions from 
Eqs. (|B5P and (|B6|) with the fits to the simulation output 
in Fig.|U Note that a^ = (3 - 2/3^) af s . This procedure 
yields a qualitatively good description of the velocity dis- 
persions produced in the simulations. Note that we do 
not assume a radial dependence of the force modification 
F (cf. Fig. [TJ to determine the velocity dispersion in 
f{R) gravity (see Fig. [3]) since the effect is blurred out 
over the different scales. In Fig. 21 we illustrate only the 
two extreme cases of the force modification, which cor- 
respond to the Newtonian and linearized f(R) scenarios, 
respectively. However, to estimate the correct halo-mass- 
dependent amplitude of the velocity dispersion, a s , in the 
transition region to the chameleon regime, we can replace 
the constant force enhancement with a weighted average 
over the modification shown in the right panel of Fig. Q] 

(cf. [H|). 
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